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Abstract 

An optimal control problem with a linear state equation and vector control with the minimization of a 
quadratic criterion is considered. Furthermore, the terminal state is fixed with fixed horizon. The relaxation 
method is proposed coupled with the shooting technique to solve this problem. We rely on shooting method 
to find initial condition of adjoint vector p(0). The convergence of the numerical procedure is presented 
where numerical experiments are also presented. We solve the problem analytically then we compare the 
results. 

Keywords: Optimal control, Shooting method, Relaxation method. 


1 Introduction 

The theory of optimal control has been well developed for over forty years. With the advances of computer 
technique, optimal control is now widely used in multi-disciplinary applications such as biological systems, 
communication networks and socio-economic systems etc. As a result, more and more people will benefit 
greatly by learning to solve the optimal control problems numerically. Realizing such growing needs, books 
on optimal control put more weight on numerical methods. In retrospect, [11] was the first and the "classic" 
book for studying the theory as well as many interesting cases (time-optimal, fuel-optimal and linear quadratic 
regulator (LQR) problems). Necessary conditions for various systems were derived and explicit solutions were 
given when possible. 

In this work, we are interested in linear quadratic optimal control problems with a fixed terminal time. 
Moreover the state is subject to some constraints and a value of the final state is fixed and vector control. 
In order to use a numerical solution, we rewrite the equations of optimality from the Pontryagin’s minimum 
principle. Then we obtain a system in which, the differential equation describing the state is equipped with 
an initial and final conditions. Note also that, the costate equation derived from the Pontryagin’s principle is 
equipped with no initial condition or terminal one to be used algorithmically. In order to determine the initial 
condition of the costate we use in this study, the shooting method [5] coupled with the relaxation method (see 
[1] and [3]). Under suitable assumptions, we analyze the convergence of the considered iterative method. 

This study is organized as follows. Section 2 defines the problem and presents the shooting method. The Sec- 
tion 3 is devoted to the description of the relaxation method coupled with the shooting one. The convergence of 
the proposed method is presented in section 4, while Section 5 contains the results of the numerical experiments. 
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2 Statement of the Problem 


Let us consider the following problem: 


Found u £ U a d such that, 
J{u) < J(v), V v £ Uad , 


Where 

T 

J[u) = l - J ((x - x d ) 2 + ku 2 )dt ; 
o 

Xd corresponding at desired state, u is the control and k control the comparative weight given to the two 
components of the cost function, i.e. , the values of k are chosen in order to grant more weight to precision or 
energy expenditure to the final state Xt fixed and minimizing the cost functional, and U a d is the closed convex 
domain. 

Under the following constraints: 

! xi = —bx i + aX 2 + Mi, 

±2 = a : Ci — bx 2 + u 2 , 
x(0) = x 0 , 
kx{T) = g, 

With G = (1, 2), 7 = 3, k= (1, 0), g = 2, T = 2, a > 0, b > 0. 

Otherwise, We can reformulate the problem as follows: 

x\ = —bx i + ax .2 + Mi, 

±2 = axi — bx 2 + M2, 

< x(0) = x 0 , 


xi(T) = 2, 


The matrix form of the state system is written as follows: 


with 



T 

/ 1 



( Xi - Xid \ 


f Ml \ 

J [(xi - x ld ,x 2 - X 2d ) 

0 

l 


x 2 - X 2 d 

+ k(ui, M 2 )* 

U2 

0 


) 


\ ) 


\ / 


If we solve the previous problem by a numerical procedure, we have then to solve an algebraic-differential 
systemOs On one hand, the state equation describing the physical system is provided with an initial condition 
Xq and the final condition Xf. On the other hand, the second equation corresponding to the costate one is 
formed with no initial condition and with no terminal condition. We will therefore use the shooting method to 
compute the value of p(0). 


2.1 Shooting method 

The shooting method is used to obtain the value of p(0) necessary to the solution of the problem characterized 
by the Pontryaguin principle. If it is possible from the condition of minimization of the Hamiltonian to express 
the extremal control with respect to (x(t),p(t)), then the extremal system is a differential system of the form 
z(t) = F(t, z(t)), where z(t) = {x(t),p(t)). With a numerical integrator beginning from zq we obtain: z*° ~ z(ti), 
where the ti are the discretized time generated by the numerical integrator. Note that in zq = (xq,Po), the value 
of xo is given by the initial condition of the problem and po denotes the value of p(0). Therefore, for different 
values of po, we will obtain the corresponding values of z z ° . We are interested in the value of z ^ ~ z(T), at the 
final time, where N is such that T = N.dt, dt being the time step, and z“$ = {x z ^,p z ^). In our case only x z £ 
are significant. Since the numerical results only depend on po they are denoted in the sequel . We define 
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the mapping G as an implicit function which gives po by using a numerical integrator such that we obtain 
x™ — Xt = 0. Let 

G : R” -> R n such that G(po) = x ^ — xt- 

Using G we define a nonlinear implicit system of n equations with n unknowns such that: 


G(po) = 0. 


(1) 


In order to solve (1), we will use the Newton’s method. The principle of the Newton’s method is described as 
follows: in the k th step, let p § be an approximation of the zero po of G ; therefore po can be written po = Po + Ap§ 
and then: 


0 = G(p 0 ) = G(pq + App) = G(pq) + 


dG 

dp 0 


(Po)-(po -Po) + o(p 0 -Po), 


which leads to the solution of the following linear system: 


dG 

dp 0 


(Po)-(Po~Po) 


-G(p k 0 ), 


where §~(Po) is the Jacobian matrix of the application po — > G(po) computed when po = p§; note that the 
mapping po — * G(po) is not explicitly known but is known numerically. So we will use a method of numerical 
derivation based on the finite difference method. To avoid the computation of §^(p k ), it is sufficient to find 
an approximation of §^(Po)- According to [4], we have two typical finite difference approximations which are 
frequently used: 


dG j 
dpoj 


(Po) 


j - 1 


t — [G j (po + hikt k ) — Gi(po + hike k )\, 


k=i 


fe= l 


or 

dG I 

-E-^iPo) « j-[Giip o + hije?) - Gj(p 0 )], 

OPOj flij 

where the hij are the given discretization step parameters of the i th equation with respect to the j th variable 
and e k are the k th vector of the canonical basis; note that, classically, we can always choose the values of h r:j 
equal to each other at a common value h. Let A ^(p, h) be a finite difference approximation, then we have: 


lim Ajj(p 0 , h) 


h^rO 


SG . . . 1 

w — [Po), i ,3 = l,...,n. 

upoj 


Let 

J (p, h) = (A^ (p 0 , h)), 

which is an approximation of the Jacobian matrix; then the approximate Newton’s method can be written as 
follows: 


Po +1 = Po - J(Po, h k )-\G(p k ), k = 0,1,..., 

The problem of convergence of this iterative process is ensured by using a result of the book of Ortega and 
Rheinbold [4]; indeed, if the discretization step hij are small enough and tend to zero then the convergence is 
ensured. 


3 Numerical Algorithm 

To solve the problem, with and without constraints, we perform the coupling of the relaxation method (see [1] 
and [3]) with the shooting method [5], the latter method being intended to calculate p(0) necessary to solve 
the differential-algebraic system obtained by applying the Pontryagin minimum principle. The steps in this 
numerical method are summarized below: 

In the unconstrained case note that x = xG ') . 

Remark 1 Steps (2), (2) and (3)correspond to the relaxation method while the following steps correspond to 
the implementation of the shooting method. 

Definition 1 An invertible matrix A is an M-matrix if A -1 > 0 and a,ij < 0 for i ^ j. 

Remark 2 M-matrices have many main properties, in particular the spectral radius of the associated Jacobi 
matrix J = I — I)~ x .A, where D is the diagonal of A, is lower than one; this property will be used in the sequel. 
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Algorithm 1 Numerical Algorithm 

0 Initial Approximation 

Control u° for t € [0,T], and of the costate variable p°(0), 
o — r <r- 0 (where r is the number of iteration), 

<0 If | u (r+1) — |> e (%where e defines the convergence threshold) do: 

- Determine the state variable x^ r \ by integration of the state equation: 

I § = Ax + Bu<- r \ 0 < t < T, 

\ 5(0) = x 0 , 

andx = Proj(x ), 


where Proj{.) is the projected operator, 

-and find the costate variable p ^ by solving: 


_ dp _ 
dt 

3 W ( 0 ), 


dH 
dx ( r ) : 


where p ^ (0) is computed by the shooting 
method, 

— Determine the control u^ r+1 ^: 


(2) 


U< r+1) i- ( U W - yS*p (r) ), 
k 


(3) 


o— Else | w.( r+1 ^ — uW | < e , 

- Determine the shooting function: 


G(p) = xW(T)-x f , 

- Solve the shooting equation by the Newton 
method and find the new value of p(0) : 

p (r+1) (0) t— p^ r \ 0) + correction , 

- r <— r + 1. 

SEnd 
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Proposition 1 If the following conditions are satisfied: 

- A is an M-matrix 

- k > ko > 0 

-p 2 (0)-p 2 (T)>0, 

then the algorithm for computing the numerical optimal control law, by the relaxation method coupled with the 
shooting method, converges whatever be the initial value of u° . 


4 Numerical Experiments 


The Hamiltonian of the system is given by: 


H(x(t),p{t),u(t),t) 


\\{xi - Xu) 2 + (x 2 - x 2 d) 2 + k{u\ + u 2 )] 

( Pl (t), P2 (t)) ( ) ( X x \ ) +{pi(t),p 2 (t)) ( H ) 

2 [(^l - Xi d) 2 + (X2 - X2d) 2 + k{u\ + v%)] +Pi(t)(-bxi + ax 2 + ui) 
p 2 (t)((ax i - bx 2 + u 2 ) +pi(t)ui +P 2 (t)u 2 - 


On the state the optimality equations can be written as follows: 

{ ^ = ^ = - h xi + ax 2 +u 1 , 

\ ^=l^= aX l~bX2 + U2, 

{ x(0) = x 0l 

< f = ^ = -bpi+ ap 2 + xi - xi d , pi(T) = ai 
\ ~^dt = = a Pl ~bp 2 +X2 -X2d, P2(T)=a 2 

oK = 0 = ku 1 +p 1 (t). 

k H =0 = ku 2 +P 2 (t). 

where on, i = 1,2 is found below. Adjunct on the previous system transversality condition: 
3 k\ ^ 0 such that: 


ip(0,x o ,T,x(T),k o ,ki) = (tfi (T, x(T)\k\), 


with 


V’l ( T ,x{T)) = kx(T) - g, 


<fi{0,xa,T,x(T),ki) = (kx(T) - g\ki). 


P(T) 


dtp 

dx{T) 


(0, xq, T, x(T), 


ki) 




(p 1 (T),p 2 (r)) = -(i,o)fc 1 = (-fc 1 ,o). 


After, the values of ai,i = 1,2 are defined, Hamilton-Pontriaguine equation’s are written as follows: 


{ x\ = —bx i + ax 2 + iii, 

±2 = ax\ - bx 2 + u 2 , 

*(0) = Xo, 

< f -pi = -bpi + ap 2 + Xi - x ld , Pi{T) = — fci 
\ ~i >2 = api - bp 2 + x 2 - x 2 d, p 2 (T) = 0 

0 = ku\ +p\{t) 

0 = ku 2 +Pi{t) 
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Then, the matrix form is given as follows: 


/ Xi \ 


/ -6 

a 

0 

0 

0 

0 

0 

\ 


( Xi 

\ 

i'2 


a 

- b 

0 

0 

1 

0 

0 



X 2 


-Pi 


-1 

0 

-6 

a 

0 

0 

0 



Pi 


-P2 

= 

0 

-1 

a 

-b 

0 

0 

0 



P2 


0 


0 

0 

1 

0 

k 

0 

0 



Ui 


0 


0 

0 

0 

1 

0 

k 

0 



U2 


v pim y 


1 0 

0 

0 

0 

0 

0 

-1 

) 


\ ki 

J 


We solve iteratively the following problem of the fixed point: 


x\ = —bx i + ax 2 + «i, 

±2 = ax i — bx 2 + u 2 , 

£i(0) = 0, £2(0) = 0, 

Pi = bpi - ap 2 -xi+ x ld , 
P 2 = -api + bp 2 - x 2 + x 2 <u 
u , = _EtW 

U1 k > 

= 

“ 2 fc > 


Pi(T) = —ki 
P2{T) = 0, 


5 Analytical Solution 

Let us derivation method at the level equation. We obtain: 

Pi = bpi - ap 2 - xi + xid, 


Pi = bpi-ap 2 -xi, 


Pi 

Pi = b(bpi - ap 2 - xi + xid) - a(—api + bp 2 - x 2 + x 2 d) - {-bx 1 + ax 2 - ~^), 
Pi = b 2 pi — abp2 — bx 1 + bxu + a 2 pi — abp 2 + ax 2 — ax 2 d + bx 1 — ax 2 + ^ , 

K 


Then, 


pi = (a 2 + b 2 + -)px ~ 2abp 2 + bxu - ax 2d , 


Pi = {a 2 + b 2 + -)pi ~ 2abp 2 + bxu - ax 2d - 


In the same way: 


P 2 = -api + bp 2 -x 2 + x 2 d 


p 2 = -api + bp 2 -x 2 , 


P2 


—a{bpi — ap 2 — X\+ x\d) + b(—api + bp 2 — x 2 + x 2d ) — {ax 1 — bx 2 



P2 


—abpi + a 2 p 2 + ax 1 — ax \ d — abpi + b 2 p 2 


bx 2 + bx 2 d 


axi + bx 2 + 


P2 

k ’ 


p 2 = -2abp! + {a 2 + b 2 + -)p 2 - axid + bx 2d , 
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Then 


P 2 = {a 2 + b 2 + -)p 2 - 2abpi - axu + bx 2 d 


( 4 ) 


derive double the equation (4), we obtain: 


Then, 


(4) leads: 


p[ 4) = (a 2 + b 2 + -)jpi - 2abp 2 , 


p^f* = (a 2 + b 2 + y)pi - 2ab(-2abp 1 + (a 2 + b 2 + y)p 2 - axid + bx 2 d), 
k k 


p i 4) = (a 2 + b 2 -)pi + 4:a 2 b 2 pi - 2ab(a 2 +b 2 + -)p 2 + 2a 2 bxid - 2ab 2 x 2 d), 


p^ = (a 2 + b 2 + j)pi + Aa 2 b 2 pi ~2cib{a 2 + b 2 + j) + 2 a 2 bxid~‘ 2 ab 2 X 2 d, 

rZ rZ 


2abp 2 = (a 2 + b 2 )pi - pi + bx i d - ax 2 d, 


By injecting (5) in (5), we obtain: 


( 5 ) 


P^ = (a 2 + b 2 + -)pi + Aa 2 b 2 pi - (a 2 + b 2 + -)((a 2 + b 2 + -)pi - pi + bxid - ax 2 d) + 2a 2 bxid - 2 ab 2 x 2 d, 


p[ 4) = (a 2 +b 2 + y)p 1 +4a 2 b 2 p 1 - ((a 2 +b 2 ) 2 + y{a 2 + b 2 ))p! + (a 2 +b 2 + ]~p! 

k k k 

— 6(a 2 + 6 2 + — )xi d H - tt(a 2 + fe 2 + — )#2 d H - 2a?bxi ( i — 2ab^X2d: 

iz iz 

p 1 4) = (2a 2 + 2 b 2 + j)pi - ((a 2 + b 2 + \) 2 - Aa 2 b 2 )p 1 + (-a 2 6 - b 3 - y + 2a 2 b)xid 

k k k 

+ (a 3 + ab 2 + y - 2 ab 2 )x 2 d, 
k 

Then, 


p[ 4) - (2a 2 + 2b 2 + ^)pi - ((a 2 + 6 2 + -j-) 2 - 4a 2 fe 2 )pi = (-b 3 - ^ + a 2 b)xu + (a 3 - afe 2 + ^)x 2 d, 

The characteristic equation corresponding in equation (6) written as: 

C 4 - 2(a 2 + b 2 + \)C 2 - (-(a 2 + b 2 + \) 2 + 4 a 2 b 2 ) = 0, 
k k 

A = 4 a 2 b 2 

The root of the characteristic equation are given by: 

1 , 


(6) 


C{ = (a 2 + b 2 + -) - 2 ab. 


C% = (a 2 + b 2 + y) + 2ab. 


Then 


Pi{t) = Xe Clt + f3e~ Clt + pe° 2t + ae~° 2t + v. 


( 7 ) 
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Determine v 1 we have: 


Pi(t) = XCie Clt — /3Cie Glt + pC^e 024 — ctC^e ° 2t , 

Pi(f) = XCfe Clt + t3Cle~ Clt + fiCle G2t + aCle~° 2 \ 

p[ 3 \t ) = A Cfe Clt - fC\e~ Glt + pCf e° 2t - aC%e~ C2 \ 

p [ 4 \t) = X Cte Clt + /3Cte~ Clt + pC%e C2t + aC$e~ C2t . 

We substitute on (6), and we obtain: 

A [C 4 - (-(a 2 + b 2 + \) 2 + 4a 2 6 2 ) - 2(a 2 + b 2 + y)C 2 ]e Glt + /3[C 4 - (-(a 2 + b 2 + \) 2 + 4a 2 6 2 ) 
k k k 

2 (a 2 + b 2 + j}C 2 ]e~ Glt + p[C 4 - (-(a 2 + b 2 + ^) 2 + 4 a 2 b 2 ) - 2(a 2 + 6 2 + \)C 2 )e G2t + a[C 4 

- (-(a 2 + b 2 + j) 2 + 4a 2 6 2 ) - 2(a 2 + b 2 + \)C'l]e~ C2t + (-(a 2 + b 2 + j) 2 + ia 2 b 2 )v 

k k k 

= {-b 3 - ^ + a 2 b)x ld + (a 3 - ab 2 + ^)x 2d , 

By identification, we obtain: 


(— b 3 
v = 

We deduce p 2 (i) of (5), we obtain: 


| + a 2 b)xid + (a 3 - ab 2 + f )x 2 d 
4 a 2 b 2 — (a 2 + b 2 + ^) 2 


P2<X> 


1 

2 ab 


[(a 2 + 6 2 + ^)pi 


Pi + kcid ^ ax 2 d2ab], 


(8) 


P 2 O) 


a 2 I l2 , 1 -t 

2a6 " [AeCl< + + ** ^ + ae ~° 2t + ^ - 2^ [AC ' eClt + P°l e ~ Glt + ^e° 2t 

aC 2 e~ G2t ) + ^~x ld - 7-rx 2d , 


P2{t) 


A [ a -±l±Ll^L ]eClt + p[ ^±^±Z^l ]e -Cit 


p[ 

( 


2 ab 

a 2 + b 2 + \- Cl C2t 


2 ab 
a 2 + b 2 + l 


2 ab 

y+i + i-ci 


2 ab 


h xu - lb*™- 


Then, 


p 2 {t) = X e Clt +/3e Clt -pe C2t -ae ° 2t + ( ° + k )i^ + -^-xi d ~ ^-x 2d - 

2 ab 2 a 2o 

Also, given the following equation: 

= bpi - ap 2 - Pi + xid, 
x 2 = — api + bp 2 - p 2 + x 2 d- 

We obtain the following results: 


(9) 


Xi(t) 


b[Xe Glt + fie Clt + pe° 2t + ae C2t + v) — afAe^ 1 * + fie Clt + pe° 2t 

a2 + b2 + K + ^x ld - ^-x 2d } - [ACie Clt - f)Cie~ Glt + pC 2 e G2t - 

2 ab 2 a 2b 


0 -C 2 t 


aC 2 e ° 2t } + x u 
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xi (t) = X(b — a — Ci)e Clt + (3(b — a + C\)e Clt + n(b + a — C 2 )e C2t + a(b + a + C 2 )e ° 2t 

r , (° 2 + b 2 + jr) 1 a 

+ [b - b ]t ,-- Xld + - X2d + Xld . 


Then, 


x\(t) = \(b — a — C\)e Clt + f3(b — a + C\)e Clt + n(b + a — C 2 )e C2t + a(b + a + C 2 )e C2t + 


b 2 -a 2 - 


2 1 


2b 


+ \x ld +^x 2d . 


( 10 ) 


x 2 (t) = —a[Xe Clt + /3e Clt + [ie C2t + ae ° 2t + v] + b[Xe Clt + fie Clt + /xe° 2< + ae ° 2t + 


+ (b 2 - a 2 - y)2 b\v + ^ xu + ^-x 2d - [ACie Clt - /3C 2 e Clt + nC 2 e C2t - aC 2 e ° 2t } + x 2d , 

k 2 Zb 


Then, 


x 2 (t) = X(b — a — Ci)e Clt + /3(b — a + C\)e Clt + /i (— b — a + C 2 )e C2t + a(— a — b — C 2 )e ° 2t 


2,1 


b 2 - a 2 + 
2 a 


-v + — x ld + -x 2d . 
Za Z 


( 11 ) 


The constants are determined by the following limits conditions: 
aii(O) = 0,x 2 (0) = 0, 

Xi(T) = 2, 

Pi(T) = -h ,p 2 (T) = 0. 

We solve the linear system numerically: 

0 = A (b — a — Ci) + fi(b — a + Ci) + fi(b + a — C 2 ) + a(b + a + C 2 ) H ^ — —v + A x\ d + -^ i x 2d , 


0 — X(b — a — Ci) + /3(b — u + Cl) + /r( — b — a + C 2 ) + a(—b — a — C 2 ) + ( 


2 6 1 2 1 2 6 
b 2 ~a 2 +i 


2 a 


)v + ^ti x ld + \x 2dj , 


2 = A (b — a — C\)e ClT + j3(b — a + C\)e ClT + n(b + a — C 2 )e C2T + a(b + a + C 2 )e C2T + 

+ \xi d + 7^X 2d , 

—k± = Xe ClT + Pe~ ClT + /ie C2T + ae~ C2T + v , 

0 = Xe ClT + Pe~ ClT - ne°*T - ae~ C2T + ± )v + _ ± X2d 

(12) is equivalent to following system: 


t.2 „2 1 


26 


b 2 -a 2 - x 


X(b — a — Ci) + f){b — a + Ci) + n(b + a — C 2 ) + a(b + a + C 2 ) - “ b k - 1 


V o'X'ld 


2b 


X(b — a — Ci) + p(b - a + C\) + n(-b - a + C 2 ) + a(-b - a - C 2 ) = -(- — - -£=xid - \x 2d , 


X (b — a — C\)e ClT + /3(b — a + C\)e ClT + /r(6 + a — C- 2 )e C2T + a(b + a + C 2 )e C2T = 2 — 

1 /y, CL 

2 Kid <2\)'^2di 

Xe ClT + /3e~ ClT + ne C2T + ae~ C2T + k\ = — u, 

Ae^ T + pe~ ClT ve c * T - ae~ C2T = - ±x ld + ±x 2d 
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Let matrix form: 


/ b — a — Ci 
b- a- Ci 
( b-a-Ci)e c ' T 

pCiT 

{ e*' 


( x \ 


/ b 2 —a 2 — -w 1 a \ 

2b V 2b X ™ 

b 2 _ Q 2 + l b 

( 2 a > V 2 2 X ™ 

p 


d 

= 

0 b 2 — a2 — j7 1 a 

^ 2b V 2 Xld 2 b X2d 

a 


—v 

V ) 


/ o, 2 ~\~b 2 -J- \ 1 1 

V ( 2 ab ) V 2 a X ^ d + 2b X ™ / 


b- 

b- 


■ a 


Ci 
a + Ci 


- a + Ci 


b 

-h 


+ CL — C2 


5.0.1 Comparison of the Two Approaches 

The numerical algorithm is implemented in Matlab. Particularly, we used the function ode 45 et / solve. In 
this example, the method converge independently with initial point ^(0) = (0.0536,0.1072). Numerical 
experience is determined by values a=l , 6 = 3 et fc = 2. Let us deduce that analytical solution and numerical 
solution correspond perfectly. For different values of k, The performance of the numerical procedure is 



Figure 1: (a) ui numerical solution (b) Vi true solution 


(a) 


(b) 



Figure 2: (a) 112 numerical solution (b) V 2 true solution 


summarized in the following table. 


k 

time C.P.U 

Iteration Number 

0.5 

1.3089 

12 

1 

2.3619 

8 

1.5 

2.7016 

7 

2 

3.7995 

7 

2.5 

5.7479 

6 


Note that the algorithm converges quickly, to a very small number of iterations required to reach convergence. 
In addition, the computing time used is very low. 
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Figure 3: (a) x\ numerical solution (b) y\ true solution 


(a) 


(b) 



Figure 4: (a) numerical solution (b) y 2 true solution 


(a) 


(b) 



Figure 5: (a) p\ numerical solution (b) q± true solution 


(a) 


(b) 



Figure 6: (a) P 2 numerical solution (b) <72 true solution 


6 Conclusions 

The goal of this study concerns the linear quadratic optimal control problems with a fixed terminal time. We 
have used the relaxation method coupled with the shooting method for solving the optimal control problem. 
The convergence of the procedure is ensured. The rate of convergence is high and the computation of time is 
fast. We concluded that relaxation method give a good results in quick time. 
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